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Abstract 

Results  are  presented  for  a  range  of  one-dimensional  shock  wave  problems  in  gaseous 
and  metallic  materials.  These  problems  were  solved  numerically  using  Flux-Corrected 
Transport  (FCT).  FCT  is  a  numerical  technique  which  achieves  high  resolution  without 
non-physical  oscillations,  especially  in  regions  of  steep  gradients  such  as  shock  fronts. 

These  types  of  problem  involve  Solving  the  Eulerian  inviscid  fluid  flow  equations,  namely 
the  continuity  equation,  conservation  of  momentum  and  conservation  of  energy,  with 
an  appropriate  equation  of  state.  For  gaseous  materials  the  ideal  gas  equation  of  state 
was  used  and  for  metallic  materials  the  “stiffened-gas”  or  the  Mie-Griineisen  equation  of 
state.  Shock  wave  problems  in  gases  included  the  one-dimensional  shock  tube  problem,  a 
shock  wave  hitting  a  density  discontinuity  and  shocks  of  equal  magnitude  colliding.  Us¬ 
ing  the  “stiffened-gas”  equation  of  state  and  the  Mie-Griineisen  equation  of  state  similar 
types  of  problems  v/ere  solved  for  metallic  materials,  for  example,  a  shock  propagating 
through  a  piece  of  metal. 

A  discussion  of  ihe  performance  of  FCT  to  accurately  model  these  problems  is  given. 
Currently  work  is  being  done  on  adding  elastic-plastic  (or  viscous)  terms  and  heat 
conduction  terms  to  the  fluid  flow  equations,  to  improve  the  description  of  flow  in  a 
solid  material.  , 
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Numerical  Modelling  of  Shocks 
in  Gases  and  Metals 


1.  Introduction 

What  happens  when  two  solids  collide?  For  impacts  which  do  not  involve  large  de¬ 
formations  a  description  of  the  behaviour  in  terms  of  elastic,  plastic  and  shock  wave 
propagation  is  needed.  The  main  aim  of  this  work  is  to  develop  a  numerical  model 
which  can  describe  the  impacts  of  solids  (involving  no  large  deformations)  in  one  and 
two  dimensions.  Before  elastic-plastic  properties  of  impacting  solids  can  be  considered 
a  thorough  understanding  of  shocks  in  solids  is  required. 

We  start  by  modelling  shocks  in  an  inviscid  gas.  Using  Flux-Corrected  Trans¬ 
port  (FCT),  the  Eulerian  inviscid  fluid  flow  equations,  governing  the  conservation  of 
mass,  momentum  and  energy  together  with  the  ideal  gas  equation  are  solved  for  a  range 
of  one-dimensional  shock  wave  problems  in  gases.  FCT  resolves  shock  fronts  over  two 
to  three  grid  points  and  moving  discontinuities  over  six  grid  points,  without  any  un¬ 
dershoots  or  overshoots.  By  choosing  a  different  equation  of  state,  which  can  describe 
the  properties  of  a  solid  material,  similar  types  of  problems  are  solved  in  a  solid.  The 
ideal  gas  equation  is  replaced  with  either  the  “stiffened-gas”  equation  of  state  or  the 
Mie-Griineisen  equation  of  state.  Treating  a  solid  as  an  inviscid  fluid  is  found  to  be 
unsatisfactory  and  the  inclusion  of  viscous  and  heat  conduction  terms  is  needed.  These 
terms  give  an  improved  description  of  shock  wave  behaviour  in  a  solid  material. 

This  report  gives  the  details  of  the  one-dimensional  problems  attempted  so  far  in 
gases  and  solids.  The  system  of  equations  and  equations  of  state  are  given  in  section  2, 
details  of  the  types  of  FCT  used  to  solve  these  equations  are  in  section  3,  a  description 
of  the  problems  attempted  is  in  section  4,  a  discussion  of  the  results  in  section  5  and 
some  concluding  remarks  are  in  section  6. 


2.  Equations 


In  this  section  the  equations  of  motion  for  a  compressible  fluid  (in  one  dimension)  are 
given,  along  with  some  simple  equations  of  state  for  gases  and  metals. 


2.1  The  system  of  equations 


In  one  dimension  the  Eulerian  inviscid  fluid  flow  equations  can  lie  written  in  conservation 
(or  vector)  form  as 


dU  dF(U) 
d  t  Ox 


(2.1) 


where 


F(U)  = 


y 


p  is  the  density,  u  the  particle  velocity,  E  the  total  energy  per  unit  volume  and  p 
the  pressure.  Heat  conduction  effects  have  been  ignored  (they  will  be  discussed  later). 
The  total  energy  per  unit  volume  can  be  expressed  as  the  sum  of  the  heat  and  kinetic 
energies, 

E  —  pi  +  \pv?  (2-2) 

where  I  is  the  internal  energy  per  unit  mass.  To  complete  this  system  of  equations  an 
appropriate  equation  of  state  needs  to  be  defined.  Usually,  the  equation  of  state  defines 
the  pressure  as  a  function  of  density  and  internal  energy, 

P  =  P{P,I)-  (2.3) 


2.2  Equations  of  state 


In  the  literature  there  is  a  wide  range  of  equations  of  state  based  on  both  theoretical  and 
experimental  results.  These  equations  can  be  extremely  complicated  and  can  describe 
a  varied  range  of  material  behaviour.  For  the  types  of  problems  solved  here  the  ideal 
gas  equation  and  simple  metal  equations  of  state  will  be  sufficient. 

For  an  ideal  gas,  equation  (2.3)  is  defined  as 

p  =  (l-l)pl  (2.4) 


where  7  is  the  ratio  of  specific  heats. 

An  equation  of  state  that  can  adequately  describe  both  gases  and  metals  is  the 
Mie-Griineisen  equation  of  state  (Ilarlow  Sc  Amsden  1971), 

P  =  PH  +  ~r»p(I  -  Ih)  (2.5) 
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where 


P.<£(1  -  *) 

[!-■*( I"  t)i2 


/h 


1/  Co(l-^)  V 


(2.6) 


(2.7) 


7,  =  2s  -  1 


(2.8) 


p#  is  the  pressure  along  the  Hugoniot  at  density  p,  ///  the  corresponding  internal 
energy  per  unit  mass,  pa  the  normal  density  of  the  material,  ca  the  speed  of  sound  in 
the  unshocked  material,  ys  the  Griineisen  ratio  and  s  the  shock  velocity  versus  particle 
velocity  slope.  The  shock  velocity  U  is  related  to  the  particle  velocity  by  the  empirical 
Hugoniot  relationship 

U  —  ca  sii.  (2.9) 

If  the  Hugoniot  cannot  be  represented  by  a  straight  line  when  plotting  shock  velocity 
versus  particle  velocity,  then  an  alternative  form  of  equation  (2.9)  has  to  be  used. 

The  expressions  for  pH  and  Ih  (equations  (2.6)  and  (2.7))  can  be  derived  from 
the  Rankine-Hugoniot  conditions  across  a  shock  front  (Hayes  1973) 

Po  j  ^  Uq 

P  u  -Uo 

P  -  Po  =  Po(U  -  u0)(u  ~  Uo) 

/  -  /0  =  }(P  +  Po)(^  -  J) 

where  the  shock  velocity  is  given  by  equation  (2.9).  The  subscript  ‘o’,  in  the  above 
equations  refers  to  the  unshocked  material. 

If  the  density  varies  only  slightly  from  the  normal  density,  then  the  Mie-Griinei- 
sen  equation  can  be  simplified  to  the  form 

P  =  clip  -  Po)  +  (7  -  1W,  (2.13) 

which  is  called  the  “stiffened-gas”  equation  of  state  (Harlow  &  Amsden  1971),  where 
7  1  =  7,.  For  large  internal  energies  this  equation  behaves  similarly  to  the  ideal  gas 

equation  (equation  (2.4)). 


(2.10) 

(2.11) 

(2.12) 
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3.  Numerical  Scheme 


A  finite  difference  approximation  to  equation  (2.1)  is  in  conservation  (or  ‘flux’)  form 
when  it  can  be  written  in  the  form, 


U 


n+1  _ 


u: 


At 

Ax 


(3.1) 


Here  U"  is  the  value  of  U  at  the  ith  grid  point  for  the  nth  time  step.  An  evenly  spaced 
grid  is  assumed.  The  Ti±x  are  called  transportive  fluxes,  and  are  functions  of  U  and 
F  at  the  nth  time  step. 


When  trying  to  model  shock  waves  numerically  we  encounter  many  difficulties. 
Some  are  due  to  the  steep  gradients  involved  (a  shock  in  an  inviscid  fluid  is  a  point 
discontinuity)  and  others  due  to  the  numerical  method  being  used,  for  example,  excessive 
implicit  diffusion  and  non-physical  oscillations.  When  choosing  a  suitable  numerical 
method  to  model  shock  waves,  we  seek  a  method  which  will  give  a  monotonic  solution 
(no  non-physical  oscillations)  and  can  accurately  model  any  discontinuities,  for  example, 
spreading  shock  fronts  (in  a  gas)  over  no  more  than  four  grid  points.  Flux-corrected 
transport  (FCT)  is  a  technique  developed  in  a  series  of  papers  by  Boris  &  Book  (1973), 
Book,  Boris  k  Hain  (1975)  and  Boris  &  Book  (1976)  which  gives  a  monotonic  solution 
with  a  high  level  of  accuracy.  To  illustrate  the  principles  of  FCT  consider  the  advection 
of  a  square  wave  (Fig.  1).  For  a  constant  velocity  field  u(x,  t)  =  u,  equation  (2.1)  reduces 
to 

H  n  fin 

(3.2) 


dp  dp 
dt+udx=0 


where  only  the  density  profile  of  the  square  wave  is  considered. 


If  a  low-order  scheme,  such  as  the  donor  cell  method  (Roache  1972),  is  used  to 
model  a  square  wave,  the  solution  (Fig.  la)  suffers  from  excessive  diffusion.  Alterna¬ 
tively,  if  a  high-order  scheme,  such  as  the  Lax-Wendroff  method  (Sod  1978),  is  used,  the 
solution  (Fig.  lb)  suffers  from  non-physical  oscillations.  When  FCT  is  used  (Figs,  lc 
and  Id,  whose  differences  of  which  will  be  discussed  later)  a  highly  accurate,  monotonic 
solution  is  obtained.  FCT  constructs  a  transportive  flux  which  is  a  weighted  average  of 
a  flux  computed  by  a  low-order  monotonic  scheme  and  a  flux  computed  by  a  high-order 
scheme.  The  weighting  procedure  blends  the  low  and  high  order  schemes  in  such  a 
way  that  the  high-order  scheme  is  used  to  as  great  as  extent  as  possible,  subject  to  the 
constraint  that  no  under  or  overshoots  are  introduced  into  the  solution. 


The  Zalesak  (1979)  FCT  procedure  (hereinafter  ZFCT)  is  as  follows: 

(1)  Compute  the  transportive  flux  given  by  some  low-order  monotonic  scheme. 

(2)  Compute  the  transportive  flux  given  by  some  high  order  (second  or  above) 

»+  j 

scheme. 
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(3)  D. 


the  ‘anti-diffusive’  fluxes: 


rp^(n) 


(4)  Compute  the  updated  low-order  (‘transported  and  diffused’)  solution: 


utd  =  u: 


A  t 

Ax 


,£(  n) 

'i+k 


,L(n) 


(5) 


Limit  the  anti-diffusive  fluxes  Ai+i  in  a  way  such  that  no  existing  extrema  are 
accentuated  and  no  new  extrema  are  created  in  the  solution  that  are  not  in  Utd  or 
Un. 


0  <  Cl+1  <  1 

—  *  '  2 


( 6 )  Apply  the  limited  antidiffusive  fluxes: 


U 


n+  1 


=  u 


td 


At 

Ax 


C(n) 

L»+  h 


C(n) 


The  critical  step  in  the  above  is  step  5,  which  is  referred  to  as  “flux-correction” 
or  “flux-limiting”.  If  the  i  coefficients  in  step  5  are  zero  then  the  low-order,  mono¬ 
tonic  solution  is  obtained,  but  if  they  are  equal  to  one  (in  regions  where  non-physical 
oscillations  are  absent)  then  the  high-order  solution  is  obtained.  In  regions  where  non¬ 
physical  oscillations  are  present  the  C!+  i  coefficients  are  calculated  in  such  a  way  that 
at  each  point  enough  implicit  nonlinear  diffusion  is  added  to  damp  out  these  oscilla¬ 
tions.  The  “flux-limiting”  process  ensures  that  U’l+1  lies  within  the  range  Um{n  to 
U max  where  Umm  and  Umar  are  determined  from  examination  of  local  Un  and  Utd 
values.  Details  of  this  ‘flux-limiting’  process  (  in  particular  how  the  C,-+i  coefficients 
are  calculated)  will  not  be  discussed  here  and  can  be  found  in  Zalesak  (1979)  FCT  can 
be  classed  as  an  implicit  artificial  viscosity  or  implicit  damping  method  since  implicit 
diffusion  or  viscosity  is  added  without  appearing  explicitly  in  the  equations  (Roache 
1972). 


Zalesak  (1979)  has  placed  the  theory  of  Boris  and  Book’s  original  FCT  (here¬ 
inafter  BBFCT)  in  a  simple,  generalised  format  which  can  accomodate  more  liberal  flux- 
limiting  techniques,  including  his  multi-dimensional  flux-limiter.  Compared  to  BBFCT, 
ZFC'T  has  an  improved  flux-limiting  stage,  but  does  not  have  the  low  phase  errors  of 
BBFCT. 

By  modifying  ZFCT  Dietachmayer  (1987)  constructed  a  scheme  that  has  the  low 
phase  errors  of  BBFCT  with  the  gererality  of  ZFCT.  This  scheme,  called  “re-ordered 
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ZFCT”(hereinafter  RZFCT),  is  based  on  ZFCT,  but  is  used  in  such  a  way  that  time 
levels  are  not  mixed  in  the  anti-diffusive  step  (see  step  6  above). 

The  procedure  is  as  follows  (Dietachmayer  1987): 

(  1  )  Compute  the  transportive  flux  given  by  some  low-order  monotonic  scheme. 

1  *  2 

(2)  Obtain  the  updated  low-order  ('transported  and  diffused’)  solution  U ' ,  using  the 
low-order  monotonic  scheme. 


u;  ~  ur 


At 

A”’ 


rj*  t- (  TI  '  rp  L(n  ) 


(3)  Compute  the  low-order  fluxes  at  (*)  level 


(  1)  Compute  the  high-order  fluxes  at  («)  level 

rH(.) 

b.  i 


(5)  Calculate  the  ‘anti-diffusive’  fluxes: 


tH(;) 


tl(V 


(6)  Limit  the  anti-diffusive  fluxes 


Ar(‘,)  C  i  A' 

7  .-  I  V  X  }  »  ~ 

I  2  *  -> 


0  '  C. 


(7)  Implement  the  antidiffusive  step. 


v?'  u; 


At 

Ar 


(  t  * )  *<■(•) 


No  extra  phase  errors  are  introduced  in  RZFC T,  making  this  procedure  more 
accurate  than  ZFCT.  Dietachmayer  (1987)  tested  a  wide  range  of  schemes,  including 
ZFCT  and  BBFCT,  and  found  RZFC  I  to  be  the  most  accurate  of  all  the  schemes. 
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Returning  to  the  square  wave  problem  discussed  earlier,  we  can  compare  the 
ZFCT  and  RZFCT  solutions  (Figs,  lc  and  Id,  respectively).  The  only  noticeable  differ¬ 
ence  between  these  solutions  is  the  smaller  amount  of  rounding  in  the  leading  edge  of 
the  square  wave  when  RZFCT  is  used  and  the  low  phase  errors  of  RZFCT  is  the  main 
reason  for  this  difference.  These  results  do  not  show  the  significant  phase  lag  of  ZFCT 
that  is  found  in  the  results  of  Dietachmayer  (1987).  This  difference  is  due  to  the  higher 
resolution  used  in  this  study.  The  ZFCT  and  RZFCT  solutions  clip  the  top  left-hand 
corner  and  bottom  right-hand  corner,  and  to  a  lesser  extent  the  top  right-hand  corner 
of  the  square  wave.  These  inaccuracies  are  due  to  the  limitations  of  the  high-order 
scheme  (Fig.  lb),  and  can  be  eliminated  by  choosing  a  more  accurate  scheme  than  the 
Lax-Wendroff  method. 


4.  Some  One-Dimensional  Test  Problems 

Some  simple  one-dimensional  problems  were  chosen  to  test  ZFCT  and  RZFCT  (see 
section  3)  for  their  accuracy  and  suitability  to  model  shock  waves. 


4.1  Gases 

Four  problems  were  used  to  test  the  ability  of  ZFCT  and  RZFCT  to  model  shock  waves 
in  a  gaseous  material  with  the  ideal  gas  equation  (equation  (2.4)).  Rusanov’s  method 
was  used  to  calculate  the  low-order  fluxes  and  the  Lax-Wendroff  method  was  used  to 
calculate  the  high-order  fluxes.  Details  of  these  numerical  methods  can  be  found  in 
Roache  (1972)  or  Sod  (1978). 

The  Shock  Tube 

The  first,  problem  was  the  one-dimensional  shock  tube  problem,  which  is  a  stringent  test 
for  the  non-linear  systems  described  by  equation  (2.1)  (see  Sod  1978).  It  is  a  simple 
problem  in  which  a  diaphragm  separates  two  regions  of  different  pressures  and  densities 
in  a  semi-infinite  tube.  The  fluid  <n  each  side  of  the  diaphragm  is  an  ideal  gas  which 
is  at  rest  at  time  t  —  0.  To  the  left  of  the  diaphragm,  the  gas  is  initially  at  a  higher 
density  and  pressure  than  on  the  right. 

Once  the  diaphragm  is  ruptured  (/  >  0),  a  rarefaction  propagates  into  the  high  - 
pressure  gas  and  a  shock  wave  followed  by  a  contact  discontinuity,  propagates  into  the 
low-pressure  gas.  Details  of  the  exact  solution  can  be  found  in  Harlow  &:  Amsden  (1971) 
and  Tyndall  (  1  988). 
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One  hundred  points  were  taken  each  side  of  the  diaphragm  (which  was  positioned 
at  x  =  0)  with  Aa'  =  0  025.  The  ratio  of  specific  heats,  7,  was  1.4  in  equation  (2.4). 
The  initial  conditions  were 

p  -  1.0  p  =  1.0  u  —  0.0  x  <  0 

and 

p  —  0.125  p  =  0.1  u  —  0.0  x  >  0. 

The  Steady  Shock 

The  second  problem  was  the  propagation  of  a  one-dimensional  shock  through  a  uniform 
polytropic  gas,  at  rest,  ahead  of  the  shock.  All  variables  downstream  of  the  shock 
were  known  and  an  initial  density  discontinuity  was  assumed.  All  the  other  variables, 
including  the  shock  velocity  U ,  were  derived  from  the  Rankine-Hugoniot  conditions  (see 
appendix  for  details).  A  schematic  diagram  of  this  problem  is  given  in  Fig.  2  and  the 
initial  conditions  are  in  Table  1. 

Shock  hitting  a  density  discontinuity 

For  the  third  problem,  the  steady  shock  described  above  was  propagated  into  a  station¬ 
ary  density  discontinuity.  This  problem  is  illustrated  in  Fig.  3  and  the  initial  conditions 
are  given  in  Table  1. 

The  collision  of  two  shocks 

This  problem  consisted  of  two  shocks  of  equal  magnitude  colliding  (see  Fig.  4).  Initial 
conditions  for  this  problem  are  given  in  Table  1. 


4.2  Metals 

The  steady  shock  problem  and  the  collision  of  two  shocks  were  the  problems  used  to  test 
the  ability  of  ZFCT  and  RZFCT  to  model  shock  waves  in  a  metallic  material.  These 
problems  involved  only  small  density  changes,  so  the  “stiffened-gas”  equation  of  state 
was  used.  Recall,  that  the  Mie-Griineisen  equation  can  be  simplified  to  the  “stiffened- 
gas”  equation  of  state  when  the  density  varies  slightly  from  the  normal  density  (see 
section  3.2).  Roe’s  upwind  scheme  (see  Roe  1981)  was  used  to  calculate  the  low-order 
fluxes  and  the  Lax-Wendroff  method  was  used  to  calculate  the  high-order  fluxes.  Three 
types  of  material  were  used  in  these  simulations,  aluminium,  copper  and  stainless  steel. 
The  physical  data  needed  for  the  Mie-Griineisen  and  “stiffened-gas”  equations  of  state 
for  each  of  these  materials  is  given  in  Table  2.  Since  they  all  behaved  similarly,  only 
the  details  and  results  for  aluminium  are  given  here. 
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The  Steady  Shock 


This  problem  is  equivalent  to  the  steady  shock  problem  in  a  gas.  A  one-dimensional  weak 
shock  was  propagated  through  aluminium  (see  Fig.  2).  All  vanables  downstream  of  the 
shock  were  specified  and  an  initial  density  discontinuity  was  assumed  (1%  compression). 
All  variables  upstream  of  the  shock  and  the  shock  velocity  U  were  derived  from  the 
Rankine-Hugoniot  conditions  in  the  same  way  as  the  steady  shock  problem  in  a  gas  (see 
appendix).  The  initial  conditions  for  this  problem  are  set  out  in  Table  3. 

The  collision  of  two  shocks 

Two  types  of  collisions  were  modelled  in  aluminium,  one  where  two  shocks  of  equal 
magnitude  collided  and  the  other  where  two  shocks  of  different  magnitudes  collided. 
The  initial  conditions  for  these  problems  are  also  set  out  in  Table  3. 


5.  Results  and  Discussion 


In  this  section  the  results  are  presented  for  the  problems  described  in  the  previous 
section. 


5.1  Gases 


The  Shock  Tube 

The  solution  to  the  shock  tube  problem  using  ZFCT  is  shown  in  Fig.  5.  It  follows 
the  general  form  of  the  exact  solution.  The  contact  discontinuity  has  been  slightly 
diffused  and  the  shock  front  has  been  well  resolved  (only  two  to  three  points  wide).  This 
solution  agrees  well  with  the  solution  given  by  Sod  (1978),  using  BBFCT  (Boris  and 
Book’s  FCT).  Neither  of  these  solutions  is  monotonic.  In  the  density  profile  there  is  an 
undershoot  at  the  end  of  the  rarefaction  (a;  =  0)  and  a  slight  overshoot  in  the  contact 
discontinuity  (f  =  1.0).  This  behaviour  is  also  present  in  the  pressure,  velocity  and 
energy  profiles.  This  non- monotonic  behaviour  is  not  caused  by  phase  errors  or  the  flux- 
correction  stage  (see  section  3)  but  by  the  initial  conditions.  If  we  integrate  numerically 
starting  from  the  exact  solution  at  some  time  t  >  0,  rather  then  rupturing  the  diaphragm 
at  t  —  0,  these  undershoots  and  overshoots  do  not  appear.  In  Fig.  6,  the  exact  solution 
at  i  0.5  was  used  as  the  initial  conditions,  then  ZFCT  was  used  to  calculate  a  solution 
at  t  -  1.0  (one  hundred  time  steps  later).  There  are  no  undershoots  at  the  end  of  the 
rarefaction  and  no  overshoots  in  the  contact  discontinuity,  and  the  shock  front  and 
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contact  discontinuity  are  the  same  widths  as  Fig.  5.  The  small  irregularities  at  x  =  0.9 
in  the  pressure  and  velocity  profiles  are  due  to  inaccuracies  in  the  numerical  scheme 
in  the  region  of  the  contact  discontinuity,  across  which  the  theoretical  pressures  and 
velocities  are  continuous. 

The  shock  tube  problem  was  also  solved  using  RZFCT,  see  Fig.  7.  Comparing 
Figs.  5  and  7,  RZFCT  gives  a  more  accurate  solution  than  ZFCT.  The  contact  discon¬ 
tinuity  is  less  rounded,  the  shock  front  is  only  two  points  wide  and  there  is  only  a  slight 
undershoot  at  the  end  of  the  rarefaction.  Using  the  one-dimensional  advection  equation, 
Dietachmayer  (1987)  tested  a  range  of  moriotonic  advection  schemes  and  found  that  the 
best  results  were  obtained  using  RZFCT.  Hence  the  remainder  of  the  problems  studied 
here  were  solved  using  RZFCT. 

The  Steady  Shock 

The  propagation  of  a  shock  through  a  gas  is  shown  in  Fig.  8.  The  shock  front  is 
represented  by  three  points  and  does  not  diffuse  as  it  travels  through  the  gas.  The 
shock  front  moved  through  the  gas  with  shock  velocity  U  =  1.19.  Apart  from  the  small 
wiggles  in  the  solution  near  x  =  0,  the  solution  is  correct  to  three  significant  figures,  but 
in  the  region  of  the  wiggles  the  solution  is  only  correct  to  twm  significant  figures.  These 
wiggles  in  the  solution  are  artifacts  of  the  initial  conditions,  which  will  be  discussed 
later. 

Shock  hitting  a  density  discontinuity 

In  Figs.  9  and  10  a  series  of  plots  show  the  interaction  of  a  shock  hitting  a  density 
discontinuity.  Shortly  before  time  t  —  1.0  (Figs.  9b  and  10b),  the  shock  (with  shock 
velocity  U  =  1.67)  collided  with  the  density  discontinuity.  If  we  look  at  the  density 
and  energy  at  t  =  1.0  (Figs.  9b  and  10b),  a  shock  has  begun  to  form  that  will  travel 
through  the  more  dense  stationary  gas  and  another  weaker  shock  has  been  reflected 
from  the  discontinuity.  The  pressure  and  particle  velocity  in  between  the  two  shocks 
should  be  constant  but,  since  the  shock  that  travels  through  the  more  dense  gas  is  still 
forming,  there  is  a  small  dip  in  the  pressure  and  the  particle  velocity  is  slightly  variable 
in  this  region.  At  later  times  t  —  1.5  and  t  =  2.5  (Figs.  9c,  9d,  10c  and  lOd)  the  shocks 
can  be  clearly  seen  and  the  pressure  and  particle  velocity  are  constant  in  this  region. 
The  shock  travelling  through  the  more  dense  gas  is  two  to  three  grid  points  wide  and 
the  reflected  shock  is  three  to  four  grid  points  wide.  Both  shocks  do  not  diffuse  as 
they  travel  through  the  gas.  The  shock  travelling  through  the  more  dense  gas  and  the 
reflected  shock  propagate  with  shock  velocities  U  —  0.97  and  U  —  —0.58,  respectively. 
The  solution  is  correct  to  three  significant  figures  except  for  the  shock  velocities  after 
the  collision,  which  are  correct  to  two  significant  figures. 

The  collision  of  two  shocks 

The  problem  of  two  colliding  shocks  is  shown  in  Fig.  11.  Before  the  shocks  collide  they 
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are  travelling  with  shock  velocities  U  —  1.67  and  U  =  —1.67.  After  the  collision  at 
x  =  0  the  shocks  are  reflected  and  travel  with  shock  velocities  U  =  1.01  and  U  =  —1.01. 
At  t  =  1.0  there  is  a  dip  region  in  the  density,  centered  at  x  =  0  and  a  slight  overshoot 
and  undershoot  in  the  particle  velocity  in  the  same  region.  These,  like  the  wiggles  in  the 
steady  shock  problem  above,  are  caused  by  the  initial  conditions,  where  the  collision 
of  the  two  shocks  can  be  considered  as  the  initial  stage  for  the  reflected  shocks.  At 
a  later  time  t  =  2.0  (Fig.  lib)  the  dip  region  in  density  transforms  into  a  region  of 
wiggles  and  the  velocity  is  monotonic.  Apart  from  these  irregularities,  the  shocks  are 
well  resolved  and  do  not  diffuse  as  they  travel  through  the  gas.  The  solution  is  correct  to 
three  significant  figures  except  for  the  irregularities  described  above,  where  the  solution 
is  correct  to  two  significant  figures. 


5.2  Metals 


The  Steady  Shock 

A  one-dimensional  weak  shock  travelling  through  aluminium, is  shown  in  Fig.  12.  The 
density  ratio  is  small  (  1%  compression)  but  the  pressure,  particle  velocity  and  energy 
ratios  are  much  larger  than  those  of  the  equivalent  gas  problem  (see  Fig.  8). Taking  this 
into  consideration,  the  shock  front  is  resolved  extremely  well,  even  though  it  is  wider 
(about  eight  points  wide)  than  the  shock  front  in  Fig.  8.  The  shock  propagates  through 
the  aluminium  with  shock  velocity  U  =  5450  ms''1,  which  is  within  0.4%  of  the  exact 
shock  velocity.  The  density,  particle  velocity  and  pressure  are  monotonic  and  involve  less 
that  0.1%  error.  However,  the  energy  is  not  monot  ic,  with  a  small  spike  at  x  =  0  m. 
When  there  is  a  larger  density  discontinuity  it  is  possible  to  see  a  corresponding  dip  in 
the  density  at  x  =  0  m.  In  the  region  of  the  spike  there  is  a  13-0%  error  in  the  energy 
and  elsewhere  there  is  less  than  0.5%  error.  In  this  problem  the  spike  does  not  affect  the 
pressure  as  the  energy  term  in  the  equation  of  state  (equation  (2.13))  is  much  smaller 
than  the  other  terms.  This  non-monotonic  behaviour  is  caused  by  the  initial  conditions, 
like  the  wiggles  and  other  irregularities  of  the  gas  problems. 

We  use  the  steady  shock  problem  to  show  how  the  initial  conditions  cause  the 
irregularities.  A  general,  linearised  solution  of  the  fluid  flow  equations  (equation  (2.1)) 
for  the  density  is 

p{x,t)  =  f(x  —  (u  +  c)t)  +  g{x  —  (u  —  c)<)  -f  h(x  —  ut)  (5.1) 

where  the  functions  /,  g  and  h  are  determined  by  the  initial  conditions.  Initially  the 
shock  was  located  at  x  =  0  m  and  assumed  to  be  an  exact  shock,  that  is,  a  point 
discontinuity  in  an  inviscid  fluid  (Fig.  13).  We  want  to  excite  a  shock  wave  which  will 
travel  with  wave  speed  lu  +  c’,  that  is,  the  functions  g  and  h  are  zero  and  /  is  non¬ 
zero  in  equation  (5.1).  When  we  model  the  theoretical  shock  numerically  there  is  a 
slight  smoothing  of  the  discontinuity  (Fig.  13)  causing  the  functions  g  and  h  to  be  small 
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rather  than  zero.  Therefore,  we  have  small  amplitude  waves  travelling  with  speeds  ‘u’ 
and  ‘u  -  c’  as  well  as  the  shock  propagating  with  speed  lu  +  c’.  The  spike  in  the  energy 
at  x  =  0  m,  and  the  wiggles  and  other  irregularities  in  the  gas  problems  are  these  small 
amplitude  waves  travelling  with  speed  lu\  Current  work  involves  including  viscous  and 
heat  conduction  terms  in  the  fluid  flow  equations,  which  should  give  a  more  accurate 
discription  of  the  initial  stage  and  diffuse  these  small  contributions.  Instead  of  assuming 
a  shock  is  a  point  discontinuity,  we  should  define  it  as  a  thin  region  where  viscosity  and 
heat  conduction  are  important. 

The  collision  of  two  shocks 

The  collision  of  two  equal  strength  shocks  in  aluminium  ,  with  shock  velocities  U  = 
5450  ms~‘  and  U  =  -5450  ms-1,  is  shown  in  Figs.  14,  15  and  16.  After  the  collision 
the  shocks  are  reflected,  similar  to  the  collision  of  two  equal  strength  shocks  in  a  gas 
(Fig.  11).  The  reflected  shocks  travel  with  shock  velocities  U  —  5430  ms-1  and  U  = 
—  5430  ms'1  which  are  within  0.4%  of  the  exact  shock  velocities.  Spikes  in  the  energy 
occur  at  the  initial  shock  positions  (r  =  — 1.5  m  and  x  —  1.5  m)  and  at  the  point 
of  collision  (x  =  0  m).  The  point  of  collision  is  where  the  shocks  change  direction, 
travelling  back  the  way  they  came,  and  can  be  thought  of  as  a  new  starting  position, 
hence  causing  a  spike  to  appear  and  also  slight  overshoots  at  the  shock  fronts.  In  the 
region  of  the  spike  at  x  —  0  m  there  is  a  6.5%  error  in  the  energy  and  in  the  region  of 
the  overshoots  a  1.2%  error. 

The  collision  of  two  shocks  of  different  strengths  in  aluminium  is  shown  in 
Figs.  17  and  18.  Before  the  collision  the  shocks  travel  with  shock  velocities  U  = 
5410  ms-1  and  U  =  -5450  ms"1  which  are  within  0.2%  and  0.4%  of  the  exact  shock 
velocities,  respectively.  After  the  collision  the  shocks  are  reflected  and  travel  with  shock 
velocities  U  =  5400  ms-1  and  U  =  —5430  ms'1,  which  are  within  0.1%  of  the  exact 
shock  velocities.  In  the  region  of  the  spikes  at  the  initial  shock  positions,  x  =  — 1.5  m 
and  x  =  2.0  m,  there  is  a  6.5%  and  13.0%  error  in  the  energy,  respectively.  At  the  point 
of  collision  x  =  0.25  m  there  is  a  2.0%  error  in  the  energy. 


6.  Conclusions 


Using  Flux-Corrected  Transport  (FCT),  the  Eulerian  inviscid  fluid  flow  equations  were 
solved  for  a  range  of  one-dimensional  shock  wave  problems  in  gases  and  metals.  Heat 
conduction  effects  were  ignored.  FCT  was  chosen  to  solve  these  problems  because  it 
gives  a  monotonic  solution  with  a  high  level  of  accuracy,  especially  in  regions  of  steep 
gradients,  such  as  shock  fronts.  FCT  resolved  the  shock  fronts  in  gases  over  no  more 
than  four  grid  points  and  contact  discontinuities  over  six  grid  points.  In  metals,  however, 
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FCT  spread  shock  fronts  over  eight  grid  points,  but,  since  the  pressure,  particle  velocity 
and  energy  ratios  were  much  larger  than  an  equivalent  gas  problem,  this  was  to  be 
expected.  In  gases  and  to  a  greater  extent  in  metals  the  numerical  solution  exhibited 
some  non-monotonic  behaviour,  found  to  be  due  to  effects  from  the  initial  conditions 
and  not  FCT.  In  all  problems  the  initial  condition  was  a  point  discontinuity,  that  is,  a 
theoretical  shock.  When  the  point  discontinuity  was  modelled  numerically  there  was  a 
slight  smoothing  of  the  discontinuity  causing  small  amplitude  waves  to  be  propagated 
through  the  solution.  The  regions  of  non-monotonic  behaviour  in  the  solutions  were 
these  small  amplitude  waves.  To  eliminate  these  small  amplitude  waves  the  initial 
condition  should  be  a  thin  region  where  viscosity  and  heat  conduction  are  important. 
Apart  from  the  regions  of  non-monotonic  behaviour,  the  numerical  solutions  for  the 
problems  in  gases  and  solids  were  within  1.0%  of  their  exact  solutions. 

Current  work  includes  adding  viscous  and  heat  conduction  terms  in  the  conser¬ 
vation  of  momentum  and  energy  equations  (equations  (2.1)),  so  that  shock  waves  in 
solids  can  be  modelled  more  accurately.  Strong  shocks  will  also  be  modelled  using  the 
Mie-Griineisen  equation  of  state  and  incorporating  elastic-plastic  terms  in  the  fluid  flow 
equations. 
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9.  Appendix 


In  this  appendix,  the  propagation  of  a  one-dimensional  shock  through  a  uniform  poly- 
tropic  gas  (see  Fig.  2)  is  considered.  Using  the  Rankine-Hugoniot  conditions  (equations 
(2.10)-(2.12)),  we  can  derive  the  variables  upstream  of  the  shock  and  the  shock  velocity. 
To  do  this,  the  variables  downstream  of  the  shock  must  be  known  and  also  one  of  the 
upstream  variables,  either  density,  pressure  or  particle  velocity. 

We  assume  an  initial  density  discontinuity  of  the  form 

Pi  =  Kp2,  (Al) 

where  K  is  an  arbitrary  constant  and  the  subscript  V  refers  to  the  unshocked  region. 
The  material  ahead  of  the  shock  is  at  rest, 


u2  =  0. 

the  Rankine-Hugoniot  conditions  become 


1  U\ 

k  =  1  ~  a 


Pl  -  P2  =  P2Uuj 


1 


h  -  I2  -  2p7 (Pi  +  P2X1  -  ft)- 

Rearranging  equation  (.A3)  gives  the  particle  velocity, 


(A2) 

(A3) 

(A4) 

(A5) 

(A6) 


From  the  ideal  gas  equation  (equation  (2.4))  the  internal  energy  per  unit  mass  can  be 
written  in  terms  of  pressure  and  density, 


(7  -  1  )p’ 


(A7) 


Substituting  equation  (A7)  into  (A5)  gives  an  expression  for  the  pressure  p\  in  terms  of 
p2,  K  and  7, 


pi_  =  +  1)  -  (7  -  1) 

P2  (7  +  1)  -  K(~i  -  1)' 

Substituting  equations  (A6)  and  (A8)  into  (A4)  gives 


(A8) 


_ 27P2-R _ 

P2H1  +  1)  -  Kfr  -  1)]’ 


(A9) 


The  pressure  of  a  very  strong  shock  rises  to  infinity  if 


Pi  _  7  +  1 

P2  7  -  1  ’ 

hence  for  the  shock  velocity  U  to  be  defined,  that  is,  U2  >  0, 

Pi  7  +  1 

P2  7  ~  1 


(410) 


(411) 


All  the  variables  upstream  of  the  shock  are  now  in  terms  of  the  initial  density 
ratio  K ,  the  pressure  in  the  unshocked  region  p2  and  the  ratio  of  specific  heats  7. 

Similar  expessions  to  equations  (-46),  (-47),  (-48)  and  (-49)  can  be  derived  for 
the  Mie-Griineisen  equation  of  state  (see  equations  (2.5)-(2.8))  and  the  “stiffened-gas” 
equation  of  state  (equation  (2.13))  when  an  initial  density  discontinuity  is  asssumed 
(equation  (-41))  and  the  variables  in  the  unshocked  region  (region  2)  are  specified. 
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Table  1 


The  initial  conditions  for  the  shock  wave  problems  in  gases,  except  the 
shock-tube  problem  (section  4.1).  The  two  regions  in  the  steady  shock 
problem  (Fig.  2)  were  la  and  2;  the  three  regions  in  the  shock  hitting 
a  density  discontinuity  (Fig.  3)  were  lb,  2  and  3;  and,  the  three  regions 
in  the  collision  of  two  shocks  (Fig.  4)  were  lb,  2  and  3b.  The  ideal  gas 
equation  (equation  (2.4))  was  used,  where  the  ratio  of  specific  heats  7 
was  1.4. 


Region  i 

Density 

Pi 

Pressure 

Pi 

Particle 

Velocity 

Ui 

la 

0.50 

0.275 

0.59 

lb 

0.25 

0.275 

0.84 

2 

0.125 

0.1 

0.0 

3a 

0.5 

0.1 

0.0 

3  b 

0.25 

0.275 

-0.84 

Table  2 


Physical  data  needed  for  the  Mie-Griineisen  and  the  “stiffened-gas”  equa¬ 
tions  of  state  (from  Matuska  &  Osborne  1985). 


Material 

Density 

p0(kgm~3) 

Griineisen 

Ratio 

7» 

Sound 

Speed 

c0(ms_1) 

Aluminium 

2710 

1.67 

5380 

Stainless  Steel 

7860 

2.46 

4610 

Copper 

8900 

1.99 

3958 
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Table  3 


The  initial  conditions  for  the  shock  wave  problems  in  aluminium  (sec 
tion  4.2).  The  two  regions  in  the  steady  shock  problem  (P4g.  2)  were 
la  and  2;  the  three  regions  in  the  collision  of  two  equal  strength  shocks 
(Fig.  4)  were  la,  2  and  3;  and,  the  three  regions  in  the  collision  of  two 
unequal  strength  shocks  (Fig.  4)  were  lb,  2  and  3. 


Region  i 

Density 

Fi(kgm~3) 

Pressure 

Pi  (Pa) 

Particle 

Velocity 

u^ras'1) 

la 

2737 

79  x  107 

54 

ib 

2724 

39  x  107 

27 

2 

2710 

0 

0 

3 

2737 

79  x  107 

-54 
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DONOR  Oclu  2  STEP  LflX-WENDRQFF 


Figure  1.  The  square-wave  test  at  time  t  =  2.0  with  A<  =  Az  =  0.01.  A  square 
wave,  which  was  initially  20  points  wide,  was  advected  at  a  constant  veloc¬ 
ity  u  —  0.2.  The  solid  line  represents  the  analytic  solution  and  the  dots 
are  the  computed  solution,  using  a)  the  donor  cell  method,  b)  the  two  step 
Lax-Wendroff  method,  c)  flux- corrected  transport  (ZFCT)  as  formulated  by 
Zalesak  (1979),  and  d)  re-ordered  Zalesak  FCT  (RZFCT)  as  formulated  by 
Dietachmayer  (198 7). 
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gure  2.  Schematic  diagram  of  the  steady  shock  problem.  A  one -dimensional  shock 
travels  at  velocity  U  through  a  uniform  gas  or  metal. 
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Figure  3.  Density  distribution  (at  time  t  =  0 )  for  a  shock  (initially  at  x  =  x0)  to  hit 
a  density  discontinuity  positioned  at  x  —  xj.  The.  density  discontinuity  is 
stationary . 
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Figure  4.  Initial  density  distibution  for  two  shocks  (of  equal  stength),  one  at  x  xa 
and  the  other  at  x  =  x,f  travelling  in  the  opposite  directions  with  absolute 
velocity  U ■ 
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Figure  5.  The  shock  tube  problem  at  time  t  ~  1.0,  with  Ax  =  0.025  and  At  =  0.005. 

The  diaphragm  was  initially  at  x  =  0.  The  circles  represent  the  solution 
given  by  ZFCT  and  the  solid  line  is  the  exact  solution.  The  profiles  for  a) 
density,  b)  pressure,  c)  particle  velocity ,  and  d)  total  energy  per  unit  volume 
are  shown. 


Figure  6.  The  shock  tube  problem  at  time  t  =  1.0,  with  Ax  =  0.025  and  A t  =  0.005. 

The  exact  solution  at  time  t  =  0.5  was  used  as  the  initial  conditions,  then 
ZFCT  was  used  to  calculate  the  solution  one  hundred  time  steps  later.  The 
profiles  for  a)  density,  b)  pressure,  c)  particle  velocity ,  and  d)  total  energy 
per  unit  volume  are  shown. 


Figure  8.  The  steady  shock  problem  (in  a  gas)  at  time  t  -  1.0,  with  At  =  0.005  and 
Ax  =  0.025.  The  RZFCT  profiles  for  a)  density,  b)  pressure,  c)  particle 
velocity,  and  d)  total  energy  per  unit  volume  are  shown.  The  shock  was 
initially  at  x  —  0. 
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Figure  9.  The  density  and  pressure  profiles  for  a  shock  hitting  a  density  discontinuity 
in  a  gas  ( using  RZFCT),  at  time  a)  t  =  0.5  and  b)  t  =  1.0.  The  density 
discontinuity  was  positioned  at  x  =  0  and  the  shock  was  initially  at  x  ~ 
-1.25. 
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Figure  11.  The  density  and  particle  velocity  profiles  for  the  collision  of  two  equal 
strength  shocks  in  a  gas  (using  RZFCT),  at  time  a)  t  —  1.0  and  b) 
t  =  2.0.  Initially  the  shocks  were  positioned  at  x  =  —1.25  and  x  =  1.25, 
with  At  =  0.005  and  Ax  =  0.025. 
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Figure  12.  The  steady  shock  problem  in  aluminium .  at  time  t  =  2.0  x  1 0 — 4  3,  with 
St  =  5.0  x  10  7  s  and  Sx  =  0.01  m.  The  RZFCT  profiles  for  a)  den- 
sity  p  (kgm  3 ),  b)  particle  velocity  u  (mf'j,  c)  pressure  p  (MPa),  and 
d)  total  energy  per  unit  volume  E  (MJm  3)  are  shown.  The  shock  was 
initially  at  x  =  0.0  m.  The  length  scale  x  is  in  metres. 
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Figure  14.  The  density  p  (kgm  3)  and  pressure  p  (MPa)  profiles  for  the  collision  of 
two  equal  strength  shocks  in  aluminium  at  time  a)  t  =  1.0  x  10~4  s  and  b) 
t  —  3.0  x  10  4  s.  The  RZFCT  solution  is  shown  where  A t  —  5.0  x  10  7  s 
and  Ax  0.01  m.  The  shocks  were  initially  at  x  =  — 1.5  m  and  x  —  1.5  m. 
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Figure  15.  Same  as  Fig.  14 ,  except  the  particle  velocity  u  (ms  1 )  and  total  energy  per 
unit,  volume  E  (MJm~3  )  at  time  a)  t  =  1.0  xl(T4  »  andb)  t  =3.0x  10  ». 
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Figure  16.  Same  as  Fig.  1J,  except  the  a)  density  p  (kgm  3 ),  b)  particle  velocity 
u  (ms  1 ),  c)  pressure  p  (MPa),  and  d)  total  energy  per  unit  volume  E 
(MJ m~3 )  at  time  i  =  4.0  x  10~4  s. 
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Figure  17.  The  density  p  (kgm~3)  and  pressure  p  (MPa)  profiles  for  the  collision  of 
two  shocks  (of  different  strengths,  see  Table  3)  in  aluminium  at  time  a) 
t  —  1.0  x  10~4  s  and  b)  t  =  3.3  x  10~4  s.  The  RZFCT  solution  is  shown 
where  A t  =  5.0  x  10-7  s  and  Ai  =  0.01  m.  The  shocks  were  initially  at 
x  =  —1.5  m  and  x  —  2.0  m. 
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colliding.  Using  the  “stiffened-gas”  equation  of  state  and  the  Mie-Griineisen  equation  of 
state  similar  types  of  problems  were  solved  for  metallic  materials,  for  example,  a  shock 
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given.  Currently  work  is  being  done  on  adding  elastic-plastic  (or  viscous)  terms  and 
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